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ABSTRACT 


Traditional  seismic  yield  estimation  is  performed  using  body  wave  magnitude  (mb)  measured  from  compressional 
wave  amplitudes  recorded  across  the  globe.  Stability  is  obtained  by  averaging  many  measurements.  These  waves 
traverse  the  earth’s  mantle,  and  are  affected  by  mantle  properties,  particularly  under  the  source  region,  as  other 
variations  are  averaged  out.  To  monitor  individual  test  sites  during  the  testing  era,  test  site  corrections  were  obtained 
by  various  means,  most  notably  the  Joint  Verification  Experiment,  and  applied  to  obtain  yield.  To  extend  yield 
estimation  to  broad  areas,  we  must  apply  an  upper  mantle  correction  on  the  fly.  We  have  investigated  two  methods 
to  map  upper  mantle  bias  over  broad  areas.  The  first  estimates  the  bias  at  individual  stations  by  inverting  for 
corrections  that  best  fit  the  collection  of  amplitudes  measured  at  all  stations.  These  measurements  were  taken  from 
global  monitoring  agencies,  including  the  International  Data  Centre  (IDC)  and  National  Earthquake  Information 
Center  (NEIC)  bulletins,  examined  separately.  We  augment  bulletin  amplitudes  by  replicating  monitoring  agency 
measurement  techniques  for  non-reporting  or  temporary  stations  in  regions  of  particular  interest,  such  as  the  Korean 
Peninsula,  including  State  University  of  New  York  (SUNY)  Binghamton's  northeast  China  and  Lamont's  Sinpo 
deployments.  The  second  method  compares  mb  to  magnitude  derived  from  regional  Lg  coda,  which  is  not  affected 
by  mantle  properties,  producing  a  map  of  the  upper  mantle  effects  across  broad  areas  where  earthquakes  occur.  The 
station-based  technique  retains  near-site  effects  that  the  event-based  technique  does  not,  thus,  resolving  any 
differences  between  the  two  techniques  is  of  great  importance. 
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OBJECTIVES 

The  North  Korea  explosion  of  October  9,  2006,  resulted  in  a  variety  of  yield  estimates  by  numerous  individuals  and 
agencies  across  the  globe.  This  situation  created  renewed  interest  in  the  technical  details  of  yield  estimation.  A 
common  means  of  estimating  yield  is  through  the  use  of  magnitude-yield  relations,  and  specifically  the  use  of 
teleseismic  P-wave  magnitude,  mb.  However,  magnitude  bias  has  long  been  recognized  as  an  issue  in  such  work, 
and  often  a  controversial  one.  Here,  we  describe  a  preliminary  look  into  magnitude  bias,  updating  some  old 
techniques  of  investigating  magnitude  bias  for  use  with  modern  data  and  databases,  with  measurements  such  as  coda 
moment  magnitudes,  and  specifically  for  the  region  of  East  Asia.  We  find  significant  potential  for  accurate 
determination  of  bias  on  a  regional  or  even  global  scale,  and  present  results  to  support  this. 

RESEARCH  ACCOMPLISHED 
Use  of  Reciprocity  to  Determine  Bias 

We  extract  a  subset  of  events  and  stations  from  the  United  States  Geologic  Survey  (USGS)-produced  monthly 
Earthquake  Data  Reports  (EDR)  catalog,  and  invert  these  data  to  determine  mb  bias  at  stations.  The  data  used  are  the 
individual  station  body  wave  magnitudes  (mbs)  reported  in  the  EDR,  during  the  period  1990-2005.  We  first 
examine  the  complete  EDR  station  magnitude  database  to  determine  which  stations  report  mb  for  a  large  number  of 
events.  We  determine  that  a  minimum  of  500  events  is  reasonable  to  preserve  good  station  statistics  and  not  greatly 
reduce  geographic  coverage.  This  results  in  a  total  of  510  stations  worldwide. 

We  then  further  reduce  the  data  set  by  restricting  the  event-station  distance  to  be  in  the  range  of  30  to  100  degrees. 
While  mb  is  often  reported  at  smaller  ranges,  the  station  measurements  at  these  small  ranges  tend  to  have  much 
greater  scatter.  We  also  eliminate  any  events  having  reported  depths  greater  than  50  km.  We  do  this  out  of  concern 
that  some  stations  may  see  a  greater  fraction  of  deep  events  relative  to  others,  and  that  this  could  create  a  spurious 
bias,  however  we  do  not  conduct  any  specific  investigation  to  see  if  this  is  really  a  problem.  We  then  limit  the  range 
of  network  mbs  to  between  4.7  and  5.7,  which  eliminates  only  about  5%  of  the  data  set.  We  use  5.7  as  an  estimate  of 
the  onset  of  mb  saturation,  and  4.7  as  an  estimate  below  which  scatter  and  imprecision  in  magnitude  increases.  Note 
that  these  limits  also  help  deal  with  data  censoring  (clipped  data  at  high  magnitudes  and  unreadable  data  at  low 
magnitudes,  problems  that  are  addressed  in  these  ranges  using  a  non-linear  maximum  likelihood  approach  to 
determine  network  magnitude,  but  which  is  not  amenable  to  linear  inversion  for  bias).  Finally,  we  limit  the  resulting 
set  of  events  to  those  that  were  recorded  at  a  minimum  of  25  stations,  so  that  particular  events  do  not  bias  individual 
stations.  The  resulting  dataset  contains  510  stations,  all  of  which  recorded  at  least  100  events;  11,609  total  events,  all 
reported  by  at  least  25  of  the  selected  stations;  and  a  total  of  647,522  station  magnitudes  (about  one  quarter  of  the 
total  number  of  EDR-reported  station  mbs). 

As  a  final  step  before  inversion,  we  eliminate  all  station  magnitude  outliers  that  are  more  than  1  magnitude  unit 
different  from  the  network  mb.  We  then  invert  the  station  mbs  using  an  LSQR  conjugate  gradient  algorithm  to 
simultaneously  determine  station  bias  and  revised  network-average  mbs,  under  the  constraint  that  mean  station  bias 
is  zero  across  the  dataset  (and  using  the  reported  network  mbs  as  the  initial  model  for  event  mb,  and  0  as  the  initial 
model  for  station  bias).  Use  of  an  inversion  technique  is  important,  since  the  station  biases  feedback  into  the 
network  mbs,  and  a  simple  mean  station  residual  could  often  be  a  poor  predictor  of  actual  station  bias  (unless,  of 
course,  those  means  were  determined  using  an  iterative  approach).  The  results  are  plotted  on  a  global  map  in  Figure 
1,  where  the  size  of  the  symbols  is  proportional  to  the  number  of  station  mbs  used  for  the  station,  and  the  color 
indicates  the  station  bias.  The  inversion  reduces  station  magnitude  residuals  by  about  15%  on  average. 
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Figure  1.  EDR  station  mb  bias. 

In  general,  expected  features  are  present.  For  example,  the  prominent  difference  in  bias  between  the  eastern  and 
western  US  is  readily  apparent.  Several  other  areas  of  thin  crust  or  active  tectonics  also  show  up  with  negative 
relative  bias.  And  ocean  island  stations  show  large  relative  positive  bias.  However,  a  significant  amount  of  scatter  is 
also  present,  as  is  evident  in  some  cases  showing  widely  differing  biases  between  nearly  co-located  stations.  In 
particular,  some  of  the  large  seismic  arrays  used  by  the  Prototype  International  Data  Center  (PIDC)  stand  out 
relative  to  nearby  stations.  It  is  known  that  the  USGS,  in  the  past,  has  used  available  amplitudes,  particularly  from 
the  PIDC,  as  mb  amplitude  measures,  but  they  were  unaware  that  these  were  not  mb  amplitudes.  This  may  have  the 
effect  of  seriously  contaminating  the  EDR  dataset  and  making  it  unsuitable  for  this  kind  of  study  without  careful 
evaluation  and  editing.  Indeed,  a  preliminary  look  at  station  residuals  as  a  function  of  time  indicates  significant  time 
dependence,  but  much  more  work  would  be  needed  to  make  such  investigations  conclusive.  Nevertheless,  the 
established  bias  at  the  Nevada  Test  Site  (NTS)  relative  to  Semipalatinsk  is  well  reproduced,  which  is  important  for 
moving  forward. 

To  avoid  the  question  of  inappropriate  or  inconsistently  determined  station  magnitudes,  we  turn  to  the  Reviewed 
Event  Bulletin  (REB)  published  by  the  PIDC  and  IDC.  These  station  magnitudes  are  determined  using  a  different 
algorithm  from  the  USGS,  but  are  very  consistent.  For  this  dataset,  we  use  the  REB  from  January  1995-May  2007. 
Given  the  relative  consistency,  and  the  smaller  number  of  stations,  we  start  with  an  initial  limit  of  100  events  per 
station  and  5  stations  per  event,  and  repeat  the  limit  of  30  to  100  degrees  in  range  and  50  km  in  depth.  We  revise  the 
mb  range  to  a  minimum  of  4.4  and  a  maximum  of  5.8,  which  is  more  appropriate  to  this  dataset,  and  eliminates 
approximately  5%  of  the  of  the  dataset,  as  the  narrower  limit  did  for  the  EDR.  This  results  in  158  stations,  15,593 
events,  and  249,506  station  mbs  for  the  inversion.  The  outlier  limit  of  1  magnitude  unit  is  again  applied,  and  then  the 
inversion  is  done  the  same  as  for  the  EDR.  The  results  are  shown  in  Figure  2. 
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Figure  2.  REB  station  mb  bias. 

The  result  for  the  REB  station  inversion  is  grossly  similar  to  that  for  the  EDR  inversion,  but  with  many  fewer  points 
(the  REB  uses  a  very  restricted  set  of  stations  relative  to  the  EDR).  However,  there  is  some  evidence  of  improved 
consistency,  particularly  in  Australia,  which  now  is  quite  uniform  and  similar  to  other  shield  areas.  This  supports  the 
inference  above  that  a  carefully  vetted  dataset  extracted  from  the  EDR  could  be  as  consistent.  Such  carefid  QC  of 
EDR-reported  station  mbs  could  also  allow  the  conservative  limits  on  numbers  of  events  and  stations  to  be  relaxed, 
greatly  increasing  the  density  of  results. 

We  are  interested  in  what  station  bias  might  reveal  regarding  magnitude  bias  in  the  vicinity  of  Korea.  We  acquired 
two  PASSCAL  datasets  in  the  vicinity  of  Korea.  One  is  the  SUNY  Binghamton  northeast  China  project  (1998-1999, 
19  stations)  and  the  other  is  Lamont’s  Sinpo  site  survey  within  North  Korea  (1995-1996,  3  stations).  The  Sinpo  data 
should  prove  quite  valuable  for  further  studies  of  North  Korea,  see  Figures  3  and  4.  The  northeast  China  dataset  is 
also  quite  interesting  and  includes  several  regional  chemical  explosions  that  were  done  as  part  of  the  experiment. 

We  have  the  full  continuous  data  from  all  22  stations.  We  select  a  number  of  events  from  global  catalogs  (EDR  and 
REB,  along  with  others)  that  would  be  at  teleseismic  range  from  these  stations  (30-100  degrees)  and  of  sufficient 
magnitude  to  be  well  recorded.  For  this  initial  survey  of  the  bias  problem,  we  analyze  data  from  three  of  the  stations: 
one  from  the  Sinpo  dataset,  HMP  (Figure  3),  and  two  from  the  northeast  China  dataset,  FROG  and  CBAI. 
Teleseismic  P  for  the  selected  events  is  then  manually  picked  at  each  of  these  stations  (Figure  4). 
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Figure  3.  Regional  earthquake  recorded  at  Sinpo  station  HMP  in  North  Korea. 
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Figure  4.  Teleseismic  earthquake  recorded  at  Sinpo  station  HMP  in  North  Korea. 

To  use  data  from  the  PASSCAL  studies  in  a  global  station  magnitude  inversion,  we  need  to  reproduce,  as  precisely 
and  accurately  as  possible,  the  computation  of  station  mb  performed  for  a  global  bulletin.  We  select  the  USGS  EDR 
primarily  due  to  its  greater  number  of  stations.  However,  reproduction  of  the  determination  of  station  mbs  for  this 
bulletin  is  not  a  trivial  exercise,  and  the  development  of  this  capability  is  important.  Using  available  documentation 
(Dewey,  et  al.,  2003,  2004)  and  a  paper  that  attempted  a  similar  reproduction  (Granville,  et  al.,  2005),  we  are  able  to 
reconstruct  the  EDR  mbs  reported  at  station  MDJ  in  eastern  China  for  a  large  number  of  teleseisms,  and  to  within 
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+/-  0.1  magnitude  units,  which  is  the  precision  to  which  station  mb  is  typically  reported  (see  Figure  5).  The  steps  in 

the  algorithm  we  use  to  reproduce  the  EDR  results  are  the  following: 

1 .  Select  a  segment  of  the  vertical  component  from  about  60  seconds  before  P  to  about  200  seconds  after; 

2.  Correct  the  segment  for  instrument  response  to  displacement  in  nm; 

3.  Filter  the  segment  using  a  1.05-2.65  Flz  2-pole  Butterworth  filter  followed  by  a  0.5-6. 5  Hz  2-pole  Butterworth 
filter; 

4.  Select  a  window  beginning  at  the  P  pick  and  extending  15  seconds  into  the  P  waveform; 

5.  Examine  the  large  window  for  the  maximum  peak-to-trough  amplitude,  using  a  sliding  window  of  1.25  seconds; 

6.  Record  half  the  maximum  peak-to-trough  amplitude  and  twice  the  time  in  seconds  between  the  positive  and 
negative  extrema  as  the  period; 

7.  Apply  a  correction  to  the  observed  amplitude  which  accounts  for  the  amplitude  attenuation  introduced  by  the 
bandpass  filters;  and 

8.  Compute  station  mb,  restricting  the  computation  to  signal-to-noise  ratio  >  2,  using  the  Gutenburg-Richter 
formula  for  mb,  along  with  the  station-event  distance  and  the  Gutenburg-Richter  attenuation  correction. 

Station  MDJ  1995-2004 


EDRM  Sta_mb  -  LANL  Sta_mb 

Figure  5.  Histogram  of  station  mb  measurements  differences  at  MDJ  between  those  reported  in  the  EDR  and 
those  computed  here. 

We  apply  this  algorithm  to  the  selected  stations  and  record  the  station  mbs.  A  cursory  examination  of  these  shows 
that  CBAI  has  a  strong  bias  for  all  results  of  nearly  1  magnitude  unit,  indicating  an  error  in  instrument  response.  We 
eliminate  the  station  rather  than  attempt  to  determine  the  source  of  the  error.  We  add  the  FROG  and  HMP  data  to  the 
EDR  dataset  used  for  Figure  1,  and  re-invert.  The  results  are  shown  in  Figure  6.  The  dataset  must  be  culled  in  a 
similar  way  to  the  original  EDR  selection,  resulting  in  a  loss  of  about  half  of  the  results. 
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Figure  6.  EDR  station  mb  bias,  with  addition  of  stations  FROG  and  HMP. 

The  remaining  station  mb  values  at  HMP  are  sufficient,  but  those  at  FROG  are  very  limited,  and  so  FROG  is 
probably  undetermined.  Interestingly,  the  results  at  HMP  indicate  that  the  station  bias  is  somewhat  shield-like,  but  it 
must  be  recognized  that  this  result  is  highly  preliminary  and  careful  work  is  needed  to  improve  the  overall  result. 
Clearly,  the  variation  shown  at  DL2  and  possibly  CN2  from  the  EDR  is  too  extreme,  with  stations  like  BJT/BJI 
showing  different  bias  while  nearly  co-located  are  probably  contaminated  by  bad  station  mbs,  etc. 


Use  of  Independent  Regional  Measures  to  Determine  Bias 

We  now  turn  to  the  use  of  regional  data  as  an  independent  estimate  of  event  size,  for  a  means  of  examining  event 
bias.  We  select  coda  moment  magnitudes  (Mw)  for  this  purpose.  We  have  a  large  number  of  network  coda  Mw 
measurement  available  for  Asia. 

We  first  investigate  the  EDR  network  mbs.  We  select  events  that  have  both  a  coda  Mw  and  a  network  mb  from  a 
global  bulletin,  subtract  the  two,  and  then  use  a  Kriging  algorithm  (Schultz,  et  ah,  1998)  to  produce  a  map  of  event 
bias.  In  this  case  we  allow  a  larger  range  of  mb  than  in  the  station  inversion  above,  permitting  a  range  of  3.8  to  5.8, 
and  we  allow  somewhat  greater  depth,  down  to  100  km.  This  is  necessary  to  keep  the  number  of  events  sufficiently 
large  for  good  kriging  statistics,  but  knowing  that  kriging  is  robust  against  outliers.  The  total  number  of  events  is 
6,417.  We  remove  the  mean,  then  decluster  the  data  to  a  grid  cell  size  of  1  degree,  and  then  krig  on  the  GRS-80 
ellipsoidal  earth,  using  an  isotropic  exponential  blending  function  with  a  range  of  10  degrees  and  a  background 
variance  of  0.4  and  a  non-bias  constraint.  The  result  is  shown  in  Figure  7.  The  fine  contours  are  0.05  magnitude 
units  of  bias,  the  heavy  contour  is  the  limit  of  resolution  at  0.2.  Small  circles  are  the  events  used.  While  the  kriging 
algorithm  ignores  continental  boundaries,  we  plot  only  the  results  within  those  boundaries,  since  the  ocean  basins 
generally  do  not  produce  accurate  coda  results.  We  restrict  the  range  of  bias  plotted  to  between  -0.5  and  +0.4  to 
maintain  a  good  range  of  color.  Those  regions  outside  this  range  are  plotted  in  magenta. 
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Figure  7.  Kriged  EDR  mb  bias  relative  to  coda  Mw. 

The  figure  shows  striking  correlations  with  crustal  and  tectonic  features.  The  thin  crust  of  China’s  southern  and 
eastern  seaboard  shows  up  as  having  relatively  positive  bias  (Mw  >  mb),  as  do  some  areas  of  recent  tectonic  activity 
such  as  the  central  tectonic  regions  of  Iran,  and  the  region  around  Bhuj  in  India.  These  regions  would  be  expected  to 
have  relatively  high  attenuation  in  the  lower  crust  and  upper  mantle,  consistent  with  relatively  low  mbs 
(more  NTS-like).  Stable  continental  regions  have  a  relatively  negative  bias  (mb  >  Mw),  such  as  the  interior  of  India 
and  north-central  Asia.  Significant  tectonic  basins  are  relatively  negative  (Sichuan,  Ordos,  and  Tarim  basins).  The 
lithosphere  tends  to  be  thick  under  these  basins,  possibly  resulting  in  lower  attenuation  in  the  lower  crust  and  upper 
mantle,  resulting  in  relatively  high  mbs  (more  Semipalatinsk-like).  The  region  around  North  Korea  is  interesting, 
showing  a  strong  gradient  from  positive  in  the  southwest  to  negative  in  the  northeast.  This  is  currently  controlled  by 
a  rather  small  number  of  events,  but  may  indicate  some  complexity  in  bias  in  this  region,  perhaps  controlled  by 
crustal  thickness  or  the  distribution  of  recent  volcanism.  The  fact  that  Punggye  would  be  in  the  region  of  relatively 
low  attenuation  (high  mb)  is  interesting  but  certainly  not  conclusive  at  this  point.  The  large  variation  in  Iran  is  also  a 
concern  for  further  investigation.  However,  the  overall  strong  correlation  with  tectonics  across  Asia  indicates  that 
this  is  a  promising  direction  for  further  study. 

Knowing  that  EDR  mbs  may  have  some  issues,  we  proceed  with  REB  mbs.  Using  the  exact  same  selection  criteria 
as  used  for  the  EDR  study,  and  the  exact  same  kriging  parameters,  we  select  6,190  events  and  present  the  result  in 
Figure  8. 
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Coda  Mw  -  REB  mb 

Figure  8.  Kriged  REB  mb  bias  relative  to  coda  Mw. 

The  resulting  pattern  is  quite  similar  to  the  EDR  result,  but  not  identical.  The  coast  of  China  appears  negative  again, 
as  does  central  Iran  and  Bhuj.  Central  India,  north-central  Asia,  the  Sichuan  basin  and  the  Tarim  basin  are  again 
negative.  The  Ordos  basin  is  almost  absent  here  relative  to  the  EDR  result,  and  the  Sichuan  basin  is  greatly 
diminished.  This  is  undoubtedly  due  to  the  subset  of  stations  used  for  the  REB,  which  may  be  having  some  bias 
effect  on  the  resulting  network  mbs.  The  gradient  in  the  vicinity  of  North  Korea  appears  to  be  a  robust  feature  of 
both  the  EDR  and  REB  mbs.  It  appears  from  the  figures  that  the  combination  of  network  averaging  of  the  station 
mbs  by  the  bulletin  authors,  and  the  kriging  is  robust  against  any  problems  in  individual  station  mbs  in  the  EDR 
bulletin. 
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CONCLUSIONS  AND  RECOMMENDATIONS 


The  results  shown  from  this  brief  exploration  of  the  mb  bias  problem  indicate  a  number  of  promising  directions  for 

further  work. 

1 .  The  station  bias  method  shows  excellent  promise  both  to  determine  broad  regional  mb  bias  and  mb  bias  at 

particular  sites  where  stations  exist.  There  are  several  aspects  of  this  method  to  pursue: 

a)  The  station  mbs  should  be  studied,  particularly  as  a  function  of  time,  to  identify  anomalous  periods  that  can 
then  be  eliminated  from  the  inversion. 

b)  AFTAC  station  mbs  should  be  used  in  an  inversion  similar  to  the  EDR  and  REB  inversions  presented. 

c)  The  algorithms  for  AFTAC’s  and  the  REB’s  station  mbs  should  be  reproduced  to  the  extent  possible  and 
applied  to  stations  not  represented  in  those  bulletins,  particularly  the  PASSCAL  datasets  used  herein. 

d)  Direct  RMS  amplitude  measurements  within  narrow  frequency  bands  and  at  a  large  number  of  global 
stations  should  be  pursued,  since  these  should  be  more  stable  than  comparatively  rough  mb  amplitude 
measures.  In  addition,  the  inversion  can  be  extended  in  this  case  to  simultaneously  operate  in  multiple 
bands  with  the  constraint  that  the  amplitude  as  a  function  of  frequency  be  linear  (that  is,  that  it  is  directly 
modeled  as  t*).  A  significant  advantage  to  such  an  inversion  is  that  t*  might  be  directly  determined  at 
events  as  well  as  at  stations. 

2.  The  Mw  bias  method  also  shows  excellent  promise  of  determining  good,  regional  maps  of  mb  bias.  Such  maps 

might  be  useful  for  predicting  bias  in  regions  without  events.  There  are  again  several  aspects  to  pursue: 

a)  Further  constraints  on  event  selection  to  improve  stability  and  resolution. 

b)  Re-computation  of  network  mbs  from  station  mbs,  but  with  strong  selection  criteria  on  the  station  mbs, 
such  as  those  discussed  above. 

c)  Use  the  map  to  predict  station  bias,  and  compare  that  to  the  results  of  the  station  bias  study.  Also,  use  it  to 
correct  station  mbs  before  re -computing  network  mbs  for  use  in  the  kriging. 

d)  Explore  the  space  of  kriging  parameters  to  see  where  the  results  are  truly  stable. 

3.  Investigate  other  independent  constraints  on  event  size,  such  as  Nuttli’s  mbLg  in  place  of  coda  Mw. 

4.  There  is  no  reason  to  restrict  the  study  to  mb,  other  magnitude  measures  (Ms,  for  example),  may  show  evidence 

of  bias  and  should  be  explored. 
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